clc;
clear;
syms x y,
% tic,
% s=int(exp(-x^2/2)*sinh(x^2+y),x,-sqrt(1-y^2),sqrt(1-y^2));
% I=int(s,y,-1,1),
% latex(I)
% vpa(I)
% toc

tic,
f=@(y,x)exp(-x.^2/2).*sinh(x.^2+y);
fh=@(y)sqrt(1-y.^2);
fl=@(y)-sqrt(1-y.^2);
I=integral2(f,-1,1,fl,fh,'RelTol',1e-20),
toc